LibrerĆ­as empleadas

library(readr)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(plotly)
## 
## Adjuntando el paquete: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

Pregunta 1: Consumo de antibióticos en España

Determinar el consumo de antibióticos por comunidades autónomas, diferenciando entre recetado y no recetado, para identificar tendencias de uso.

Pregunta 2: Relación entre PIB y resistencia a antibióticos en la UE

Analizar si existe una relación significativa entre el Producto Interno Bruto de los países de la Unión Europea y los niveles de resistencia a antibióticos.

pibPP <- read_table("INPUT/DATA/datos_pib.tsv")
## Warning: Missing column names filled in: 'X14' [14]
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   `freq,na_item,ppp_cat,geo\TIME_PERIOD` = col_character(),
##   `2012` = col_double(),
##   `2013` = col_double(),
##   `2014` = col_double(),
##   `2015` = col_double(),
##   `2016` = col_double(),
##   `2017` = col_double(),
##   `2018` = col_double(),
##   `2019` = col_double(),
##   `2020` = col_double(),
##   `2021` = col_double(),
##   `2022` = col_double(),
##   `2023` = col_double(),
##   X14 = col_logical()
## )
View(pibPP)

str(pibPP)
## spc_tbl_ [42 Ɨ 14] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ freq,na_item,ppp_cat,geo\TIME_PERIOD: chr [1:42] "A,VI_PPS_EU27_2020_HAB,GDP,AL" "A,VI_PPS_EU27_2020_HAB,GDP,AT" "A,VI_PPS_EU27_2020_HAB,GDP,BA" "A,VI_PPS_EU27_2020_HAB,GDP,BE" ...
##  $ 2012                                 : num [1:42] 30 133 30 121 47 170 91 84 124 128 ...
##  $ 2013                                 : num [1:42] 29 133 31 121 46 171 84 86 125 130 ...
##  $ 2014                                 : num [1:42] 30 132 30 121 47 171 81 88 127 129 ...
##  $ 2015                                 : num [1:42] 30 131 31 121 48 171 83 89 124 128 ...
##  $ 2016                                 : num [1:42] 30 130 31 120 49 166 88 89 125 128 ...
##  $ 2017                                 : num [1:42] 30 127 31 118 50 160 90 91 124 130 ...
##  $ 2018                                 : num [1:42] 30 128 32 118 52 158 91 92 124 129 ...
##  $ 2019                                 : num [1:42] 30 126 33 118 53 153 93 93 121 126 ...
##  $ 2020                                 : num [1:42] 31 125 33 118 55 155 91 93 123 133 ...
##  $ 2021                                 : num [1:42] 31 122 34 120 57 157 94 92 120 135 ...
##  $ 2022                                 : num [1:42] 34 124 35 120 62 159 94 90 117 136 ...
##  $ 2023                                 : num [1:42] 35 123 35 118 64 154 95 91 115 127 ...
##  $ X14                                  : logi [1:42] NA NA NA NA NA NA ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   `freq,na_item,ppp_cat,geo\TIME_PERIOD` = col_character(),
##   ..   `2012` = col_double(),
##   ..   `2013` = col_double(),
##   ..   `2014` = col_double(),
##   ..   `2015` = col_double(),
##   ..   `2016` = col_double(),
##   ..   `2017` = col_double(),
##   ..   `2018` = col_double(),
##   ..   `2019` = col_double(),
##   ..   `2020` = col_double(),
##   ..   `2021` = col_double(),
##   ..   `2022` = col_double(),
##   ..   `2023` = col_double(),
##   ..   X14 = col_logical()
##   .. )
# Cambio de nombre columna
colnames(pibPP)[1] <- "pais"

# Nos quedamos con las Ćŗltimas letras
pibPP$pais <- substr(pibPP$pais, nchar(pibPP$pais) - 1, nchar(pibPP$pais))
View(pibPP)

lista_pais <- list("BE", "BG", "CZ", "DK", "DE", "EE", "IE", "EL", "ES", "FR", "HR", "IT", "CY", "LV", "LT", "LU", "HU", "MT", "NL",
               "AT", "PL", "PT", "RO", "SI", "SK", "FI", "SE")

# nos quedamos solo con los paĆ­ses de la UE


pib <- pibPP %>% filter(pais %in% unlist(lista_pais))
View(pib)


# quitar la columna nula

pib <- pib[, colSums(is.na(pib)) < nrow(pib)]
pib
## # A tibble: 27 Ɨ 13
##    pais  `2012` `2013` `2014` `2015` `2016` `2017` `2018` `2019` `2020` `2021`
##    <chr>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>
##  1 AT       133    133    132    131    130    127    128    126    125    122
##  2 BE       121    121    121    121    120    118    118    118    118    120
##  3 BG        47     46     47     48     49     50     52     53     55     57
##  4 CY        91     84     81     83     88     90     91     93     91     94
##  5 CZ        84     86     88     89     89     91     92     93     93     92
##  6 DE       124    125    127    124    125    124    124    121    123    120
##  7 DK       128    130    129    128    128    130    129    126    133    135
##  8 EE        74     76     78     76     77     79     82     83     85     86
##  9 EL        71     72     72     70     68     67     66     66     62     63
## 10 ES        91     90     90     91     92     93     91     91     83     84
## # ℹ 17 more rows
## # ℹ 2 more variables: `2022` <dbl>, `2023` <dbl>
pib_2022 <- pib %>% select(pais, `2022`)
View(pib_2022)


# ggplot del pib en el 2022 (primero lo ponemos en descendente)

pib_2022_desc <- pib_2022 %>% arrange(desc(`2022`))


# sustituir las etiquetas de los paĆ­ses

pib_2022_desc <- pib_2022_desc %>%
  mutate(pais = case_when(
    pais == "SK" ~ "Eslovaquia",
    pais == "SI" ~ "Slovenia",
    pais == "EE" ~ "Estonia",
    pais == "MT" ~ "Malta",
    pais == "LV" ~ "Latvia",
    pais == "HR" ~ "Croatia",
    pais == "EL" ~ "Greece",
    pais == "BG" ~ "Bulgaria",
    TRUE ~ pais # Mantiene los nombres que no estƔn en la lista
  ))

Tras el tratamiento de estos datos se puede ver de manera grƔfica la comparativa entre los paƭses.

grafico_pib <- ggplot(pib_2022_desc, aes(x = reorder(pais, -`2022`), y = `2022`)) +
  geom_bar(stat = "identity", fill = "gold") +
  labs(x = "PaĆ­s", y = "Valor en 2022", title = "PIB por PaĆ­s en 2022") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # Rotar etiquetas para mejor legibilidad

grafico_pib

Una vez observada la primera variable, generamos los datos del valor medio de la resistencia a antimicrobianos en los países europeos para comprobar si existe una relación o, efectivamente no tiene que ver.

Incidencia_enfermedades <- read_delim("INPUT/DATA/ECDC_encuesta_AMR_incidencia_enfermedades.csv", 
                      delim = ",", escape_double = FALSE, trim_ws = TRUE)
## Rows: 98334 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (8): HealthTopic, Population, Distribution, Unit, RegionCode, RegionName...
## dbl (2): Time, CategoryIndex
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Incidencia_enfermedades
## # A tibble: 98,334 Ɨ 10
##    HealthTopic         Population Distribution Unit   Time RegionCode RegionName
##    <chr>               <chr>      <chr>        <chr> <dbl> <chr>      <chr>     
##  1 Antimicrobial resi… Acinetoba… R - resista… %      2012 AT         Austria   
##  2 Antimicrobial resi… Acinetoba… R - resista… %      2012 AT         Austria   
##  3 Antimicrobial resi… Acinetoba… R - resista… %      2012 AT         Austria   
##  4 Antimicrobial resi… Acinetoba… R - resista… %      2012 AT         Austria   
##  5 Antimicrobial resi… Acinetoba… R - resista… %      2012 BE         Belgium   
##  6 Antimicrobial resi… Acinetoba… R - resista… %      2012 BE         Belgium   
##  7 Antimicrobial resi… Acinetoba… R - resista… %      2012 BE         Belgium   
##  8 Antimicrobial resi… Acinetoba… R - resista… %      2012 BE         Belgium   
##  9 Antimicrobial resi… Acinetoba… R - resista… %      2012 BG         Bulgaria  
## 10 Antimicrobial resi… Acinetoba… R - resista… %      2012 BG         Bulgaria  
## # ℹ 98,324 more rows
## # ℹ 3 more variables: CategoryIndex <dbl>, Category <chr>, Value <chr>
summary(Incidencia_enfermedades)
##  HealthTopic         Population        Distribution           Unit          
##  Length:98334       Length:98334       Length:98334       Length:98334      
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##       Time       RegionCode         RegionName        CategoryIndex  
##  Min.   :2000   Length:98334       Length:98334       Min.   :1.000  
##  1st Qu.:2008   Class :character   Class :character   1st Qu.:1.000  
##  Median :2013   Mode  :character   Mode  :character   Median :2.000  
##  Mean   :2013                                         Mean   :2.158  
##  3rd Qu.:2018                                         3rd Qu.:3.000  
##  Max.   :2022                                         Max.   :4.000  
##    Category            Value          
##  Length:98334       Length:98334      
##  Class :character   Class :character  
##  Mode  :character   Mode  :character  
##                                       
##                                       
## 
View(Incidencia_enfermedades) 

incidencia_2022 <- Incidencia_enfermedades %>%
  filter(Time == 2022) 
  
incidencia_2022MF <- incidencia_2022 %>%  
  filter(Category == 'Male' | Category == 'Female')

View(incidencia_2022MF)


incidencia_2022MF$Value <- as.numeric(incidencia_2022MF$Value)
## Warning: NAs introducidos por coerción
media_poblacion_genero <- incidencia_2022MF %>%
  arrange(RegionCode, Population) %>%  # Ordenar los datos
  group_by(RegionCode, Population) %>%
  mutate(mean_value = (Value + lead(Value)) / 2) %>%
  ungroup()

incidencia <- media_poblacion_genero %>%
  mutate(grupo = substr(Population, 1, 3)) %>%
  select(-Unit, -Category, -CategoryIndex, -Value, -Distribution)%>%
  filter(!is.na(mean_value))

# Ver el dataframe
print(incidencia)
## # A tibble: 754 Ɨ 7
##    HealthTopic           Population  Time RegionCode RegionName mean_value grupo
##    <chr>                 <chr>      <dbl> <chr>      <chr>           <dbl> <chr>
##  1 Antimicrobial resist… Acinetoba…  2022 AT         Austria         0     Aci  
##  2 Antimicrobial resist… Acinetoba…  2022 AT         Austria         3.10  Aci  
##  3 Antimicrobial resist… Acinetoba…  2022 AT         Austria         0     Aci  
##  4 Antimicrobial resist… Acinetoba…  2022 AT         Austria         0.909 Aci  
##  5 Antimicrobial resist… Enterococ…  2022 AT         Austria         3.58  Ent  
##  6 Antimicrobial resist… Enterococ…  2022 AT         Austria         8.53  Ent  
##  7 Antimicrobial resist… Enterococ…  2022 AT         Austria         0     Ent  
##  8 Antimicrobial resist… Enterococ…  2022 AT         Austria        90.1   Ent  
##  9 Antimicrobial resist… Enterococ…  2022 AT         Austria        30.8   Ent  
## 10 Antimicrobial resist… Enterococ…  2022 AT         Austria         2.49  Ent  
## # ℹ 744 more rows
View(incidencia)


media_por_bacteria <- incidencia %>%
  mutate(grupo = substr(Population, 1, 3)) %>%                # Crear una columna con las primeras letras
  arrange(grupo, RegionCode) %>%                              # Ordenar por grupo y RegionCode
  group_by(grupo, RegionCode) %>%                             # Agrupar por grupo y RegionCode
  summarise(media_mean_value = mean(mean_value, na.rm = TRUE)) %>%  # Calcular la media para cada combinación de grupo y RegionCode
  ungroup()
## `summarise()` has grouped output by 'grupo'. You can override using the
## `.groups` argument.
View(media_por_bacteria)

# Dataframe final con el valor de media final de todas las bacterias
media_por_pais <- media_por_bacteria %>%
  group_by(RegionCode) %>%
  summarise(media = mean(media_mean_value, na.rm = TRUE))

Gracias a estos datos podemos estudiar las bacterias que generan resistencia de manera individual, y ver su incidencia en los paĆ­ses.

# grƔfico que te dice quƩ bacteria tiene mƔs incidencia
ggplot(media_por_bacteria, aes(x = grupo, y = media_mean_value)) +
  geom_boxplot() +
  labs(title = "Distribution of Mean Incidence by Bacteria Group",
       x = "Bacteria Group", y = "Mean Incidence") +
  theme_minimal()

# grƔfico que te dice quƩ paƭses tienen esa media para cada bacteria

# Crear el conjunto de datos interactivo con highlight_key
incidencia_keyed <- highlight_key(media_por_bacteria, ~RegionCode)

# Crear el grƔfico ggplot con el texto configurado para el tooltip
scatter_plot <- ggplot(incidencia_keyed, aes(x = grupo, y = media_mean_value, color = RegionCode, text = paste("PaĆ­s:", RegionCode))) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Population vs Mean Incidence Value by Region",
       x = "Population", y = "Mean Incidence Value") +
  theme_classic()

# Convertir el grƔfico ggplot en un grƔfico interactivo de plotly con highlight
interactive_scatter_plot <- ggplotly(scatter_plot, tooltip = "text") %>%
  highlight(on = "plotly_click", off = "plotly_doubleclick", color = "red", opacityDim = 0.2)
## `geom_smooth()` using formula = 'y ~ x'
# cuando haces click en un paƭs te dice solo los puntitos de ese paƭs, cuando haces dobleckick en otra parte del grƔfico, se desaparece.

# Mostrar el grƔfico interactivo
interactive_scatter_plot

Pregunta 3: Consumo y resistencia en sectores especĆ­ficos en la UE

Se pretende explorar cómo el consumo de antibióticos en diferentes sectores se relaciona con los niveles de resistencia a antibióticos en diferentes países de la UE.

Funcion para sacar los datos JSON Tenemos la carpeta ā€œResistencia_Antibioticos_UEā€ con 27 archivos .jsonld por lo que hemos creado 3 funciones para obtener los datos. Estos datos se corresponden a la ganaderĆ­a de la UE.

La funcion obtener_nombre obtiene una lista con los nombres de los archivos que contiene la carpeta

La función obtener_archivo carga cada archivo, accede al enlace necesario y descarga la carpeta.zip (con los datos) de cada archivo y las guarda en una carpeta

La función leer_archivo recorre la carpeta que contiene los archivos.zip los descomprime y carga los datos excel. Hace un dataframe con los datos y los almacena en el global enviroment de cada uno de los paises.

library(httr)
library(tidyverse)
library(jsonlite)
library(readxl)
library(rjson)

obtener_nombre<-function(carpeta){
  archivos <- as.list(list.files(path =carpeta))
  lista_nombres<-list()
  for(i in 1:length(archivos)){
    posicion1<-regexpr("_", archivos[[i]])
    posicion2<-regexpr("\\.", archivos[[i]])
    subcadena<-substr(archivos[[i]], posicion1+1, posicion2-1)
    lista_nombres[[i]]<-subcadena
  }
  
  return(lista_nombres)
  
}


obtener_archivo<-function(direccion){
  lista_paises<-obtener_nombre(direccion)
  lista_enlace<-list()
  direccion_archivos<-list()
  
  for(i in lista_paises){
    cada_pais<-paste0("AMR_",i,".json")
    lista_enlace[i]<-cada_pais
  }
  for(i in lista_enlace){
    
    cada_archivo<-paste0(direccion,"/",i)
    direccion_archivos[i]<-cada_archivo
  }
  
  
  for(i in direccion_archivos){
    pais<-fromJSON(file= i)
    enlace<-pais$links$archive
    respuesta_archivo <- GET(enlace)# Hacer la solicitud HTTP para descargar el archivo
    nombre_archivo<-basename(enlace)#Extrae el nombre del archivo de la URL
    
    if (status_code(respuesta_archivo) == 200) {# Verificar si la solicitud fue exitosa (código 200, código estÔndar HTTP que significa "OK")
      # Guardar el archivo ZIP localmente en formato binario
      writeBin(content(respuesta_archivo, "raw"), nombre_archivo)
      print("Archivo ZIP descargado correctamente.")
      unzip(nombre_archivo, exdir = "carpeta_destino", overwrite = TRUE)
    } else {
      print(paste("Error al descargar el archivo. Código de respuesta:", status_code(respuesta_archivo)))
    }
    
   
  }
  
}
obtener_archivo("INPUT/DATA/Resistecia_Antibioticos_UE")






leer_archivo <- function(carpeta) {
  carpeta_destino <- carpeta
  archivos_zip <- list.files(carpeta_destino, pattern = "\\.zip$", full.names = TRUE)
  
  # Iterar sobre los archivos .zip y descomprimirlos
  for (archivo in archivos_zip) {
    # Descomprimir el archivo .zip
    archivos_extraidos <- unzip(archivo, exdir = carpeta_destino, overwrite = TRUE)
    print(paste("Descomprimido:", archivo))  # Imprimir cada archivo que se descomprime
    
    # Filtrar el archivo .xlsx entre los extraĆ­dos
    archivo_xlsx <- archivos_extraidos[grepl("\\.xlsx$", archivos_extraidos)]#AquĆ­, grepl() busca archivos cuyos nombres terminen con .xlsx 
    
    # Verificar si hay algĆŗn archivo .xlsx descomprimido
    if (length(archivo_xlsx) > 0) {
      # Leer el archivo Excel como dataframe
      datos_xlsx <- read_excel(archivo_xlsx[1])  # Leer el primer archivo .xlsx encontrado
      
      # Asignar el dataframe al Global Environment usando el nombre del archivo como variable
      nombre_variable <- make.names(basename(archivo_xlsx[1]))  # Crear un nombre de variable vƔlido
      assign(nombre_variable, datos_xlsx, envir = .GlobalEnv)  # Asignar el dataframe al Global Environment
      
    } 
  }
}


leer_archivo("carpeta_destino")

Juntar los dataframes que tenemos en el enviroment 1. Listar los nombres de todos los dataframes que terminan en ā€œ.xlsxā€ (ajusta si es necesario) 2. Convertir los nombres a una lista de dataframes usando mget() 3. Unir todos los dataframes en uno solo usando bind_rows

eval=FALSE

nombres_dataframes <- ls(pattern = "_AMR_PUB\\.xlsx$")


lista_dataframes <- mget(nombres_dataframes)

 
df_combinado <- bind_rows(lista_dataframes, .id = "origen")

Cargamos el global enviroment En el que tenemos guardado todos los dataframe de cada pais resultante de la funcion leer_archivo, y el dataframe df_combinado del codigo anterior

load("OUTPUT/paisesUE_datos_AMR.RData")

Modificación y obtencion del dataframe con los datos seleccionamos las columnas que vamos a necesitar y filtramos en la columna de las bacterias solo las patogénicas.

Cambiamos los nombres de las columnas y lo asignamos al dataframe

paises_UE_df<-df_combinado%>%
  select(rep_Country_name,rep_Country_code,zoonosis_name,matrix_name,totUnitsTested,totUnitsPositive,sampUnitType_name,sampType_name,MIC_name,CUTOFFVALUE)%>%
  mutate(zoonosis_name = sub(".*", "", zoonosis_name))%>% # Extraer solo la primera palabra
  filter(zoonosis_name != "Escherichia coli, non-pathogenic, unspecified")


nuevos_nombres <- c("NombrePais", "Codigo", "zoonosis_name","OrigenMuestra", "TotalMuestras","MuestraPositiva","Tipo_Unidad_Muestra","TipoMuestra","MIC_name","ValorCorte")  # Modifica segĆŗn el nĆŗmero de columnas

colnames(paises_UE_df) <- nuevos_nombres

Mapa Interactivo de positivos en ganaderia por paises de la UE

paises_UE <- c(
  "Cyprus", "France", "Lithuania", "Czechia", "Germany", 
  "Estonia", "Latvia", "Sweden", "Finland", "Luxembourg", 
  "Belgium", "Spain", "Denmark", "Romania", "Hungary", 
  "Slovakia", "Poland", "Ireland", "Greece", "Austria", 
  "Italy", "Netherlands", "Croatia", "Slovenia", "Bulgaria", 
  "Portugal", "Malta"
)



mapa_mudo <- st_read("INPUT/DATA/mapaMundi")  
## Reading layer `ne_10m_admin_0_countries' from data source 
##   `C:\Users\deyan\fuentes\seminario_fuentes\INPUT\DATA\mapaMundi' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 258 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.6341
## Geodetic CRS:  WGS 84
mapa_mundo_europa <- mapa_mudo %>% 
  filter(NAME %in% paises_UE)

mapa_mundo_europa <- mapa_mudo %>% filter(NAME %in% paises_UE)


positivos_por_ciudad <- paises_UE_df %>%
  group_by(NombrePais) %>%
  summarize(total_pruebas = sum(TotalMuestras, na.rm = TRUE),
            total_positivos= sum(MuestraPositiva, na.rm = TRUE)) %>%
  mutate(ratio_positivo = (total_positivos / total_pruebas) * 100)

# Unir datos de positividad al mapa
mapa_mundo_europa$NAME <- as.character(mapa_mundo_europa$NAME)
positivos_por_ciudad$NombrePais <- as.character(positivos_por_ciudad$NombrePais)

# Realizar el join usando las columnas correctas
mapa_mundo_europa <- mapa_mundo_europa %>% 
  left_join(positivos_por_ciudad, by = c("NAME" = "NombrePais"))



mapa <- ggplot(mapa_mundo_europa) +
  geom_sf(aes(fill = ratio_positivo)) +
  scale_fill_viridis_c(option = "plasma", na.value = "gray") +
  labs(title = "Tasa de Positividad por PaĆ­s en Europa",
       fill = "Tasa de Positividad (%)") +
  coord_sf(xlim = c(-30, 50), ylim = c(35, 72), expand = FALSE) +  # Ajustar lĆ­mites para hacer zoom en Europa
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))


# Convertir el grƔfico a interactivo con plotly
mapa_interactivo <- ggplotly(mapa)

mapa_interactivo